Ground state properties of the 2D disordered Hubbard model 
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We study the ground state of the two-dimensional (2D) disordered Hubbard model by means of 
the projector quantum Monte Carlo (PQMC) method. This approach allows us to investigate the 
ground state properties of this model for lattice sizes up to 10 x 10, at quarter filling, for a broad 
range of interaction and disorder strengths. Our results show that the ground state of this system 
of spin-1/2 fermions remains localised in the presence of the short-ranged Hubbard interaction. 
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I. INTRODUCTION 

The electronic transport properties of disordered sys- 
tems have been the subject of much investigation in 
physics. From the pioneering work of Anderson (1958) 0] 
it is known that in three dimensions (3D) the eigenstates 
of a non-interacting electron gas in a random potential 
become localised at the Fermi energy above a critical 
value of the disorder strength Wc- In this regime, the 
eigenstates decay exponentially in space and hence can- 
not carry a current; thus the system is an insulator. For 
disorder strength, W, lower than Wc, the eigenstates are 
extended and diffusive transport takes place in the sys- 
tem in accordance with Ohm's law. However, for two- 
dimensional (2D) systems, it was shown by the scaling 
theory of Abrahams et al |^ that all states are localised 
for any disorder strength. Thus, it appears that there is 
no metal-insulator transition (MIT) for non-interacting 
electrons in 2D. The properties of non- interacting elec- 
trons in random potentials have since been studied sys- 
tematically and the main physical effects have been un- 
derstood 1^. In this context, the experimental observa- 
tion by Kravchenko et a/ Q of a transition from insulat- 
ing to metallic behaviour, as seen in the resistivity as a 
function of temperature of the 2D electron gas, came as 
a great surprise to the community. The existence of this 
transition from insulating to metallic behaviour as a func- 
tion of density has been confirmed by other groups ■ 
The experiments were carried out on very high mobility, 
low electron density {us) samples which correspond to a 
regime where the electron-electron interactions {E^e) are 
much stronger than the Fermi energy (Ep), such that the 
dimensionless parameter rs(w Eee/Ep) lies in the range 
5-50. This indicates the importance of electron-electron 
interactions in these systems. Further, it was shown ex- 
perimentally that the application of an in-plane magnetic 
field (Bp) drives the system insulating Q. Since such a 
field can only couple to the spin, this experiment indi- 
cates the important role played by the spin degrees of 
freedom. While the early experiments have stimulated 
a spate of new experimental results, there has been no 
satisfactory theoretical explanation of the phenomenon 



of the 2D MIT, to date. 

The effects of interactions in disordered systems have 
been studied from the metallic side in great detail, where 
the interactions are relatively weak |l^. The study of 
interactions in the localised phase were mainly carried 
out within the mean-field approximation, which led to 
a number of important results, for example, the Efros- 
Shklovskii gap in the density of states near the Fermi 
level PI , but could not take into account quantum inter- 
ference effects, important in this many-body system. The 
investigation of a simple model of two interacting parti- 
cles (TIP) in the localised phase showed that short-range 
attractive/repulsive interactions can lead to destruction 
of localisation and propagation of pairs of particles on 
a length scale much larger than the one-particle localisa- 
tion length ||l^. Thus, the effects of interaction on the lo- 
calised phase are non-trivial and deserve a detailed study. 
This, however, is not an easy task. Indeed, even the an- 
alytical expressions for the matrix elements of the inter- 
action in the localised phase are not known hence, 
numerical studies of the problem become important. 

Recent numerical approaches to the question have in- 
cluded the studies of persistent currents by exact diag- 
onalisation of small 2D clusters ||lj] and Hartree-Fock 
based calculations without [ p^ and with residual interac- 
tion [ p^JT?! . These approaches led to some interesting in- 
dications but did not allow the study of sufficiently large 
systems and/or sufficiently many particles. Other ap- 
proaches based on level spacing statistics of many-body 
states made possible the study of larger systems and 
showed the existence of ergodic (delocalised) states for 
low energy excitations, but not at the ground state [ p^ . 
All these studies were carried out for spinless fermions. 
Recently, the properties of fermions with spin on a disor- 
dered 2D lattice were investigated using a finite temper- 
ature quantum Monte Carlo (QMC) method |19|. The 
temperature dependence of the resistivity obtained nu- 
merically indicated a transition from insulating to metal- 
lic behaviour for sufficiently strong interactions and weak 
disorder strength. However, these calculations were car- 
ried out at finite temperatures and technical problems 
(" fcrmion sign problem" ) did not permit the analysis 
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of the ground state. This is not completely relevant to 
the experiments which were carried out at temperatures 
much below the Fermi energy and thus require a bet- 
ter understanding of the properties of the ground state, 
as in the general scenario of quantum phase transitions. 

To investigate the properties of the ground state of 
the disordered interacting fermionic system, we choose 
the Hubbard model with site diagonal disorder. This 
could be considered as an important first step on the 
way to investigations of more complicated models with 
Coulomb interactions, which might be more appropriate 
for experiments at low densities. We study the ground 
state of our model on the square lattice by the projector 
QMC (PQMC) method. We use different characteristics 
to investigate the extent of the ground state wavefunction 
for a broad range of model parameters: disorder strength 
(W), interaction strength ([/) and filling factor (z/). The 
studies were carried out in the Sz ~ sector, with equal 
numbers of particles with up and down spins. This is 
thus the first numerical study of the ground state of a 
disordered, interacting system of fermions with spin. 

This paper is organised as follows. In the next sec- 
tion, we describe the model and the method used. In 
the third section we present our results for the averaged 
Green function, the charge densities and the inverse par- 
ticipation ratios, all of which we use to characterise the 
ground state. We present a summary of our conclusions 
in the last section. 

II. MODEL AND METHOD 

The two-dimensional disordered Hubbard model on a 
square lattice is given by 

H = Ha + Hi 

' i,t7 ' ^ i 

where the ('li.o-) are the creation (annihilation) op- 
erators for a fermion of spin a at site i with periodic 
boundary conditions, "hi^ is the number operator for spin 
a at site z, t is the hopping parameter, the Hubbard pa- 
rameter, measures the strength of the screened inter- 
action and Ci, the energy of site j is a random number 
drawn from a uniform distribution [—W/2,W/2], which 
parametrizes the disorder. The first two terms represent 
the Anderson Hamiltonian and the last term represents 
the interaction Hj. In the limit W = 0, this Hamiltonian 
reduces to the usual Hubbard model. The filling factor 
ly = Np/{2 X 7V^), where Np is the number of fermions 
(particles) and A'' the linear dimension of the system; 
thus the total number of sites is iV^ . 

We obtain the ground state properties of this model 
by the PQMC method. The PQMC method was initially 
developed for the Hubbard model and has been used to 



obtain reliable results for large lattices . The method 
can be generalised in a straightforward manner to include 
random site energies. We now present some details of the 
calculation for completeness and refer the reader to the 
literature for more detailed accounts. 

The PQMC method consists in obtaining the true 
ground state |'0o) of the Hamiltonian(l) by projection 
from a trial wavefunction |0) that is not orthogonal to 
the true ground state of the system, 

IV'o)^lim " 1'^^ . (2) 
^(^|e-2eH|0) 

The trial wavefunction is usually formed from the eigen- 
states of the non-interacting Hamiltonian (orbitals filled 
up to the Fermi level). In this case, we choose the eigen- 
states of the Hamiltonian Ha, thus including the random 
potential. To carry out Monte Carlo (MC) simulations of 
this quantum Hamiltonian, it is first necessary to map it 
onto an effective classical Hamiltonian. Thus, the projec- 
tion operator exp{—QH) is first Trotter decomposed as 

(exp{~ATHA)exp{~ATHi)^ with = At x L. This 

introduces a systematic error of order (Ar)^ due to non- 
commutation of Ha and Hj. The interaction is then de- 
coupled by a discrete Hubbard-Stratonovich (H-S) trans- 
formation, by the introduction of N'^ x L Ising-like fields. 
Since the complete summation over these degrees of free- 
dom is too time-consuming to be practical, the method 
reduces to a MC sampling of physical properties, which is 
the second source of error, the statistical error. It is im- 
portant to note that during the MC process, each config- 
uration of Ising spins is assigned a weight, which is inter- 
preted as a probability. This quantity is positive definite 
only at half-filling for the uniform Hubbard model. The 
problems that arise from the non-positive-definite nature 
of this quantity are referred to as the 'fermion sign prob- 
lem' in the literature and are known to be particularly 
severe slightly away from half-filling in finite temperature 
methods and restrict the lowest temperature that can be 
attained in the simulation (in the clean limit). 

We have studied system sizes of up to 10 x 10 at quarter 
and one-eighth fillings (50 and 25 particles). We carried 
out extensive checks on the MC parameters to assure 
ourselves of convergence, as described in Ref. |2^, but 
in the presence of disorder. Thus, we chose Q — 3.0, 
with L = 30, to have At = 0.1. This, with the symmet- 
ric Trotter decomposition reduces the systematic error 
to (At)3 « 0.001. We checked for statistical convergence 
of our data in several ways. By varying the number of 
sweeps after equilibration, we determined that 1000 MC 
sweeps are sufficient for equilibration and 2000 further 
sweeps for property estimates. We carried out measure- 
ments until the standard deviations on our values were 
of the order of the systematic error. We also tested our 
results against results obtained from exact diagonalisa- 
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tions for small system sizes and the results for the charge 
density, rii = riii + riii = '^^{nia-)^ presented in Fig. 
1 show good agreement with the exact results. Conver- 
gence is of course the best for the ground state energy 
as compared to other physical quantities, and we have a 
relative accuracy of 10~^ when compared to exact calcu- 
lations. As for the effect of disorder on the sign problem, 
it was possible to study the 10 x 10 lattice for U/t = 6 for 
disorder strengths W/t of up to 7-10. Our measure of the 
severity of the sign problem is to consider the quantity f 
= 1- (number of negative determinants)/ (total number 
of determinants). In all the cases considered, we have f 
= 0.999, which indicates that the sign problem is under 
control. 

2 I ■ 1 ■ 1 ■ 1 ■ 1 ■ 1 



1.5 - 




FIG. 1. Comparison of exact diagonalisation (circles) and 
PQMC (squares) results for the charge density (n^) per site 
(i) for a 2-chain Hubbard model, dotted lines correspond to 
the upper chain and dashed lines to the lower chain, with 
system size 6x2, U/t = 2, W/t = 10 and filling at 4 particles. 

From the simulations, it is possible to obtain ground 
state expectation values of the single particle Green func- 
tion, Gij = CT^-o")' where the average is a MC av- 
erage. Further, we can obtain ground state expectation 
values of other one- and two-body operators, such as the 
charge density and the charge-charge correlation func- 
tions. Each disorder realisation constitutes a full PQMC 
calculation. The properties are averaged over 16 disor- 
der realisations. Thus, we have obtained the evolution 
of the Green function with distance, the charge densities 
and the inverse participation ratios. Our results will be 
described in the following section. 

III. RESULTS AND DISCUSSION 

To characterise the properties of the ground state, we 
study the correlation function defined as 



'^W = ^2 <E K^T^i.T) + (4i«i,i>l'^i-.>), (3) 

where r = i — j is the vector in the plane between the sites 
labelled i and j, and the averages are carried out over 
the ground state eigenfunction and the different disorder 
realisations for all possible initial positions of r, i.e. all 
corresponding i and j. C(r) is simply related to the 
Green function dj already defined in Section II. With 
this definition, C(0) = iV"^(E.,("*T + "u)^) ^ 4;^^ in 
the limit of weak disorder and C(0) ^ 4;/ in the strongly 
localised limit. The dependence of C(r) on distance r is 
related to the localisation of the eigenstate, i.e. we expect 
exponential decay of this quantity for localised states and 
slow decay at long distances for extended states. We also 
study the direction averaged correlation function C{r) 
which now depends only on the distance r = |r|. 
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FIG. 2. The decay of the direction averaged correlation 
function C{r) vs. r for a 30 x 30 system for (7/t = 0, with W/t 
= 2 (circles), 7 (squares), 10 (diamonds), 15 (up-triangles) 
and a 50 X 50 system at W/t — 15 (down-triangles), averaged 
over 16 disorder realisations, at quarter filling. 

In Fig. 2 we show the decay of C{r) with r, for the 2D 
Anderson model, system (1) nt U/t = 0, for various dis- 
order strengths. The change from flat behaviour with r 
at weak disorder, when the eigcnfunctions are delocalised 
in the finite sized system (W/t = 2), to asymptotic expo- 
nential decay for stronger disorder, W/t > 10, (when the 
localisation length is smaller than the system size) , is ev- 
ident. We note that the initial non-exponential decay in 
the localised case is due to the fact that the ground state 
eigenfunction is a superposition of one-particle eigen- 
states of different energies. In fact, the one-particle lo- 
calisation length li of a state depends on its energy and 
therefore the many-body state, which is the Slater de- 
terminant of the 1-particle states up to the Fermi level, 
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initially decays more rapidly. This is due to the low- 
lying states of smaller localisation lengths and it is only 
in the asymptotic limit that the decay of the many-body 
ground-state is determined by the maximum li at Ep. 
This physical structure of many-body states complicates 
the observation of asymptotic exponential decay corre- 
sponding to the one-particle localisation length Ii{Ef) at 
the Fermi level. Despite these complications, the asymp- 
totic slope is seen to depend strongly on W (Fig. 2), 
which is consistent with the exponential growth of li with 
decreasing W in 2D In view of this, the investigation 
of the correlation function C{r) in the presence of inter- 
actions should tell us the impact of interactions on the 
localisation properties of eigenstates. 
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FIG. 3. Decay of C(r) for Hamiltonian(l) with 
U/t = 0, W/t = 7,iV = 10 and = 1/4 (50 fermions on a 
10 X 10 lattice), averaged over 16 W values. The upper part(a) 
shows the decay in 3D form for the interval < C(r) < 0.01, 
the lower part(b) is a contour plot of the same data. 



The dependence of C(r) on r, shown in Figs. 3a,b 
for the 2D disordered, non-interacting Hubbard model 
model [U/t = 0), also clearly shows localisation of the 
ground state eigenfunction. The decay (as seen from the 
contour plot Fig. 3b) is approximately symmetric in r 
which is due to averaging over different disorder realiza- 
tions. Hence, it should be useful to study this quantity 
also in the interacting case, to clarify the ground state 
properties. 
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FIG. 4. Same as in Fig. 3a,b with U/t = 6 and the same 
disorder realisations. 

The behaviour of C(r) as a function of r, for rela- 
tively strong interaction strength, {U = 6t) is shown in 
Figs. 4a, b. The comparison with the non-interacting case 
(Figs. 3a, b) clearly shows that even such a strong inter- 
action produces only a slight change in the correlation 
function. We observe similar behaviour for other disorder 
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and interaction strengths ( 5 < W/t < 10, < U/t < 6, 
1/ = 1/4, 1/8, results not presented here). This, in our 
opinion, provides direct evidence that even in the pres- 
ence of strong interactions, the ground state of the system 
remains locahsed. This conclusion is futher supported 
by the data for the direction averaged correlation func- 
tion, C{r) presented in Fig. 5. In fact, this direction 
average further smoothens fluctuations due to disorder. 
Indeed, even the introduction of relatively strong inter- 
actions {U/t ~ 6) affects this function very weakly. 
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FIG. 5. The decay of the direction averaged correlation 
function C{r) vs. r for the Hamiltonian (1) with U/t — 
(circles), 2 (squares), 6 (diamonds) for W/t = 7,N = 10 and 
filling = 1/4 (50 fermions on a 10 x 10 lattice), averaged 
over the same 16 disorder realisations. 

As an alternative test for localisation, we use another, 
more indirect method, similar to the approach presented 
in Ref. As discussed in we vary the ampli- 

tude of on-site disorder Ci for sites i along one vertical 
line of the square lattice, as — > a x with a = 1.1 
and 1.3, corresponding to 10% and 30% change in dis- 
order. We then study the charge density difference Sp^ 
produced by this perturbation, as a function of distance 
X from the original line. We average over all sites with 
the same x and additionally average log \Spx \ over 16 dis- 
order realisations. The comparison of data for 10% and 
30% variation shows that we are in the linear response 
regime. The results are presented in Fig. 6. For U/t — 0, 
the response function shows a sharp drop from the ini- 
tial peak followed by a slower decay at longer distances. 
This behaviour is qualitatively similar to the decay of 
the correlation function (Fig. 2, Fig. 5), for the same 
physical reasons as analysed above. Introducing interac- 
tions doesn't at all affect the main structure of the curve 
which drops very quickly from the centre by more than 



one order of magnitude. We interpret this as a sign of a 
localised ground state. At the same time, we note that 
there is a slight difference introduced by interactions at 
the tails of the response functions. However this corre- 
sponds to a density variation less than 0.1%, which is at 
the limit of the accuracy of our calculation. In the light 
of the ensemble of data, we conclude that the ground 
state in the presence of interactions remains localised. 

-1.5 I ' : ' : ' : ' : ' 1 
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FIG. 6. Behaviour of (log \5p^\) with x for U/t = (filled 
symbols), 2 (open symbols) and W/t = 2 (squares) and 7 
(circles) with A'' = 10, v = 1/4 and a = 1.1, averaged over 
the same 16 realisations. 

In a sense the Hubbard repulsion leads to local rear- 
rangements of charge and does not seem to influence the 
long-distance properties of the system. This is clearly 
illustrated in Figs. 7a,b where the interaction leads to a 
more homogeneous charge distribution (note the change 
in the vertical scale) but doesn't drastically change the 
global profile. 

This point of view is further borne out by the data for 
the inverse participation ratio (IPR), presented in Fig. 
8. The IPR, deflned as ^ = (E« E.("0'))> 
gives the average number of sites visited per particle. ^ 
is bounded from above at fixed filling, with ^ < S^max — 
(2j/)^^, corresponding to the weak disorder limit and 
from below with ^ > = 0.5 in the strongly localised 
limit. We note that C(0) x ^ = 2i^. We have studied the 
IPR as a function of system size and interaction strength 
at 1/4 (Fig. 8) and 1/8 (Fig. 9) fiUing. At U/t = 0, the 
IPR naturally increases with decreasing disorder as states 
become more extended. At strong disorder {W/t > 5) 
^ is not sensitive to the system size since the states are 
localised and at fixed filling the IPR is counted per parti- 
cle. The Hubbard interaction smoothly increases the IPR 
but does not introduce a qualitative change. This corre- 
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spends to the local reorganisation of charge introduced 
by U leading to a more homogeneous charge density dis- 
tribution, as discussed above. The data presented for 
1/8 filling in Fig. 9 show qualitatively similar behaviour. 
However, the size variation is more restricted in this case 
( < 10 in our studies). 




FIG. 7. Charge density (rii) at site i for U/t — (upper 
figure (a)) and 6 (lower figure (b)) for a 10 x 10 lattice, W/t 
= 7, z.=l/4. 
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FIG. 8. IPR (^) vs. U/t for system sizes 10 x 10 (circles), 
8x8 (squares) and 6x6 (diamonds) and increasing disor- 
der strength from top to bottom, W/t — 0.5 (dot-dashed), 
2 (long-dashed), 5 (dashed), 7 (dotted) and 10 (solid) lines, 
at quarter filling, averaged over the same 16 disorder realisa- 
tions. 
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FIG. 9. Same as in Fig. 8 for for system sizes 8x8 (circles) 
and 4x4 (squares) at 1/8 filling. 

At this point, it is interesting to compare our studies 
with a recent finite temperature QMC study of a similar 
model Q . These authors considered the Hubbard model 
on a square lattice, with off-diagonal disorder, in contrast 
to our study. This is due to the fact that the method 
used, a finite temperature QMC method, suffers from a 
severe sign problem in the presence of diagonal disorder. 
For off-diagonal disorder, the situation becomes better. 
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but the problem persists, restricting the lowest accessible 
temperatures. Their studies of several physical charac- 
teristics, including the conductivity (obtained by approx- 
imate analytic continuation of the imaginary time Green 
function), indicate the presence of an interaction induced 
metal-insulator transition in their model. However, this 
method is not adapted to analysis of the ground state 
properties. Our results are not in direct contradiction 
with this study, since it is fully possible that the ground 
state remains localised, while the low-lying excited states 
become delocalised. Such a situation has been observed 
in numerical studies of spinless fermions with Coulomb 
interactions on a 2D lattice with disorder . Our result 
directly demonstrates the localised nature of the ground 
state even in the presence of strong interactions. This 
is important in the framework of general studies of zero- 
temperature quantum phase transitions. 



IV. CONCLUSIONS 

We have used the projector quantum Monte Carlo 
(PQMC) method to study the ground state of the 2D 
disordered Hubbard model. This method allows us to 
study systems of up to 50 spin-1/2 fermions on a 10 x 10 
lattice, for interaction strengths U /t up to 6 and a broad 
interval of disorder strengths. The comparison of several 
properties in the absence and presence of the Hubbard 
interaction allows us to conclude that interactions lead 
to local rearrangements of charge but do not destroy the 
localised structure of the ground state, within the range 
of parameter values studied. 

These results indicate that short-range interactions are 
probably insufficient to bring about a quantum phase 
transition in the ground state of this system. Thus, it 
becomes important to consider the effect of long-range 
Coulomb interactions for electrons on a disordered lat- 
tice. 

We thank M-B. Lepetit for useful discussions and the 
IDRIS, Orsay, for allocation of resources on their super- 
computers. 
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